A constrained scheme for Einstein equations based on Dirac gauge and spherical 

coordinates 
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We propose a new formulation for 3+1 numerical relativity, based on a constrained scheme and 
a generalization of Dirac gauge to spherical coordinates. This is made possible thanks to the 
introduction of a flat 3-metric on the spatial hypersurfaces t — const, which corresponds to the 
asymptotic structure of the physical 3-metric induced by the spacetime metric. Thanks to the joint 
use of Dirac gauge, maximal slicing and spherical components of tensor fields, the ten Einstein 
equations are reduced to a system of five quasi-linear elliptic equations (including the Hamiltonian 
and momentum constraints) coupled to two quasi-linear scalar wave equations. The remaining three 
degrees of freedom are fixed by the Dirac gauge. Indeed this gauge allows a direct computation of 
the spherical components of the conformal metric from the two scalar potentials which obey the 
wave equations. We present some numerical evolution of 3-D gravitational wave spacetimes which 
demonstrates the stability of the proposed scheme. 
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I. INTRODUCTION AND MOTIVATIONS 

Motivated by the construction of the detectors LIGO, 
GEO600, TAMA and VIRGO, as well as by the space 
project LISA, numerical studies of gravitational wave 
sources are numerous (see [1, 2] for recent reviews). The 
majority of them are performed within the framework 
of the so-called 3+1 formalism of general relativity, also 
called Cauchy formulation, in which the spacetime is foli- 
ated by a family of spacelike hypersurfaces. We propose 
here a new strategy within this formalism, based on a 
constrained scheme and spherical coordinates, which is 
motivated as follows. 



A. Motivations for a constrained scheme 

In the 3+1 formalism, the Einstein equations are de- 
composed in a set of four constraint equations and a set 
of six dynamical equations [1, 3]. The constraint equa- 
tions give rise to elliptic (or sometime parabolic) partial 
differential equations (PDE), whereas the PDE type of 
the dynamical equations depends on the choice of the co- 
ordinate system. Various strategies can then be contem- 
plated: (i) free evolution scheme: solving the constraint 
equations only to get the initial data and performing the 
time evolution via the dynamical equations, without en- 
forcing the constraints; (ii) partially constrained scheme: 



using some of the constraints to compute some of the 
metric components during the evolution and (iii) fully 
constrained scheme: solving the four constraint equations 
at each time step. 

In the eighties, partially constrained schemes, with 
only the Hamiltonian constraint enforced, have been 
widely used in 2-D (axisymmetric) computations (e.g. 
Bardeen and Piran [4], Stark and Piran [5], Evans [6]). 
Still in the 2-D axisymmetric case, fully constrained 
schemes have been used by Evans [7] and Shapiro and 
Teukolsky [8] for non-rotating spacetimes, and by Abra- 
hams, Cook, Shapiro and Teukolsky [9] for rotating ones. 
We also notice that the recent (2+l)+l axisymmetric 
code of Choptuik et al. [10] is based on a constrained 
scheme too. 

Regarding the 3-D case, almost all numerical studies 
to date are based on free evolution schemes^. It turned 
out that the free evolution scheme directly applied to 
the standard 3+1 equations (sometimes called ADM for- 
mulation) failed due to the development of constraint- 
violating modes. An impressive amounts of works have 
then been devoted these last years to finding stable evolu- 
tion schemes (see [12] for an extensive review and [13] for 
a very recent work in this area). Among them, a large 
number of authors have tried to introduce coordinates 
and auxiliary variables so that the dynamical equations 
become a first-order symmetric hyperbolic system. How- 
ever these approaches have revealed very limited success 
in practice. Another approach has become very popular 
in the last few years: the so-called BSSN formulation, 
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originally devised by Shibata and Nakamura [14] and re- 
introduced by Baumgarte and Shapiro [15]. It has shown 
a much improved stability with respect to the standard 
ADM formulation. Indeed the most successful compu- 
tations in numerical relativity to date are based on that 
formulation (e.g. [16, 17]). 

All the approaches mentioned above favor first-order 
hyperbolic equations with respect to elliptic equations. 
In particular, they employ a free-evolution scheme, avoid- 
ing to solve the (elliptic) constraint equations. The main 
reason is neither mathematical nor physical, but rather 
a technical one: for most numerical techniques, solving 
elliptic equations is CPU time expensive. In this article, 
we present an approach which is based on the opposite 
strategy, namely to use as much as possible elliptic equa- 
tions and as few hyperbolic equations as possible. More 
precisely we propose to use a fully constrained-evolution 
scheme and to solve the minimum number of hyperbolic 
equations: the two wave equations corresponding to the 
two degrees of freedom of the gravitational field. The 
main advantages of this procedure are that (i) elliptic 
equations are much more stable than hyperbolic ones, 
in particular their mathematical well-posedness is usu- 
ally established, (ii) the constraint-violating modes that 
plague the free-evolution schemes do not exist by con- 
struction in a fully constrained evolution, (iii) the equa- 
tions describing stationary spacetimes are usually elliptic 
and are naturally recovered when taking the steady-state 
limit of the proposed scheme. Besides, let us point that 
some very efficient (i.e. requiring a modest CPU time) 
numerical techniques (based on spectral methods) are 
now available to solve elliptic equations [18, 19]. Very 
recently some scheme has been proposed in which the 
constraints, re-written as time evolution equations, are 
satisfied up to the time discretization error [20]. On the 
contrary, our scheme guarantees that the constraints are 
fulfilled within the precision of the space discretization 
error (which can have a much better accuracy, thanks to 
the use of spectral methods). 

To achieve this aim, we use maximal slicing, as long 
as a generalization of Dirac gauge to curvilinear coordi- 
nates. This gaTige fixes the spatial coordinates (x') in 
each hypersurface t ~ const. It has been introduced by 
Dirac in 1959 [21] as a way to fix the coordinates in the 
Hamiltonian formulation of general relativity, prior to 
its quantization (see [22] for a discussion). Dirac gauge 
has been discussed in the context of numerical relativity 
first by Smarr and York, in their search for a radiation 
gauge in general relativity [23]. But they disregarded it 
as being not covariant under coordinate transformation 
(a;*) I— »• (a;* ) in the hypersurface t = const. They pre- 
ferred the minimal distortion gauge, which is fully covari- 
ant and allows for an arbitrary choice of the coordinates 
(x*) in the initial hypersurface. Here we show that if one 
introduces a flat 3-metric on each spatial hypersurface, in 
addition to the physical 3-metric induced by the space- 
time metric, the Dirac gauge can be made covariant. This 
enables the use of curvilinear coordinates, whereas Dirac 



original formulation was only for Cartesian coordinates. 
However, contrary to the minimal distortion gauge, this 
generalized Dirac gauge still determines fully the coor- 
dinates in the initial slice (up to some inner boundary 
conditions if the slice contains some holes). 



B. Motivations for spherical coordinates 

Since the astrophysical objects we want to model (neu- 
tron stars and black holes) have spherical topology, it is 
natural to use spherical coordinates (x*) = {r,9,(p) to 
describe them. In particular, spherical coordinates and 
spherical components of tensor fields enable one to treat 
properly the boundary conditions (i) at the surface of 
fluid stars, as well as at some black hole (apparent) hori- 
zon, and (ii) at spatial infinity or at the edge of the com- 
putational domain. For a binary system, two systems of 
spherical coordinates (each centered on one of the ob- 
jects) have proved to be successful in the treatment of 
binary neutron stars [24] and binary black holes [25]. 



1. Outer boundary conditions 

For elliptic equations, spherical coordinates allow a 
natural 1/r compactification which permits to impose 
boundary conditions at spatial infinity [19, 26]. In this 
way, the imposed boundary conditions are exact. 

For wave equations from a central source, a spherical 
boundary of the numerical domain of integration allows 
to set non-reflecting boundary conditions [27]. Moreover 
the use of spherical components of the metric tensor al- 
lows, in the Dirac gauge, an easy extraction of the wave 
components. This results from the asymptotic transverse 
and traceless (TT) behavior of Dirac gauge and the fact 
that a TT tensor wave propagating in the radial direction 
is well described with spherical components. 



2. Black hole excision 

Spherical coordinates clearly facilitate black hole exci- 
sion. Moreover for stationary problems, one has usually 
to set the lapse to zero on some sphere r — const, in or- 
der to preserve the time-independent behavior of slicing 
of stationary spacetimes [28, 29]. As we discuss in Ap- 
pendix A, using spherical components of the metric ten- 
sor and shift vector is crucial is setting boundary condi- 
tion on an excised 2-spherc with vanishing lapse function. 
In fact, because of the degeneracy of the operator acting 
on the above quantities when the lapse is zero, one can 
impose boundary conditions on certain components, and 
not on the others. In Cartesian components (i.e. linear 
combinations of spherical components) , the imposition of 
boundary conditions could not be done simply. 
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3. Fulfilling the Dirac gauge 

We will show that, when using spherical coordinates, 

the Dirac gauge condition can be imposed easily on spher- 
ical components of the metric tensor. Indeed, we propose 
to use the Dirac gauge to compute directly some met- 
ric components from the other ones. This seems difficult 
with Cartesian components (even with spherical coordi- 
nates). 

4- Spherical coordinates and numerical techniques 

Despite the above strong advantages and although they 
have been widely used for 2-D (axisymmetric) compu- 
tations [4-9, 30-32], spherical coordinates are not well 
spread in 3-D numerical relativity. A few exceptions are 
the time evolution of pure gravitational wave spacetimes 
by Nakamura et al. [30] ^ and the attempts of comput- 
ing 3-D stellar core collapse by Stark [33]. This situation 
is mostly due to the massive \isagc of finite difference 
methods, which have difficulties to treat the coordinate 
singularities on the axis ^ = and 6 = n, and at the 
origin r = 0. On the contrary, spectral methods em- 
ployed mostly in our group [19, 34] and Cornell group 
[35] , deal without any difficulty with the singularities in- 
herent to spherical coordinates. Let us note that in other 
fields of numerical simulation, like stellar hydrodynam- 
ics, spherical coordinates are well spread, for instance in 
the treatment of supernovae [36, 37] . 

C. Plan of the paper 

Wc; start the present study by introducing in Sec. II a 
conformal decomposition of the 3-1-1 Einstein equations 
which is fully covariant with respect to a background 
flat metric. This differs slightly from previous conformal 
decompositions (e.g. [14, 15]) by the fact that our con- 
formal metric is a genuine tensor field, and not a tensor 
density. Then in Sec. Ill we re-write the conformal 3-1-1 
Einstein equations in terms of the covariant derivative 
with respect to the; flat background metric. This enables 
us to introduce the (generalized) Dirac gauge in Sec. IV 
and to simplify accordingly the equations. We introduce 
as the basic object of our formulation the difference h be- 
tween the inverse conformal metric and the inverse flat 
metric. At the end of Sect. IV, we present an explicit 
wave equation for h. In Sec. V, we introduce spherical 
coordinates and explicit the equations in terms of ten- 
sor components with respect to an orthonormal spherical 
frame. We show how the Dirac gauge can then be used 
to deduce some metric components from the others in 



a quasi-algebraic way. The resolution of the dynamical 
3-|-l equations is then reduced to the resolution of two 
(scalar) wave equations. A numerical application is pre- 
sented in Sec. VI, where it is shown that the proposed 
scheme can evolve stably pure gravitational wave space- 
times. Finally Sec. VII gives the concluding remarks. 
This article is intended to be followed by another study 
which focuses on the treatment of boundary conditions 
at black hole horizon(s). Here we present only in Ap- 
pendix A a preliminary discussion about the type and 
the number of inner boundary conditions for black hole 
spacetimes. 



II. COVARIANT 3-1-1 CONFORMAL 
DECOMPOSITION 

A. 3+1 formalism 

We refer the reader to [1] and [3] for an introduction to 
the 3+1 formalism of general relativity. Here we simply 
summarize a few key equations, in order mainly to fix the 
notations^. The spacetime (or at least the part of it under 
study) is foliated by a family of spacelike hypersurfaces 
St, labeled by the time coordinate t. We denote by n 
the future directed unit normal to Ej. By definition n, 
considered as a 1-form, is parallel to the gradient of t: 

n = -Ndt. (1) 

The proportionality factor N is called the lapse function. 
It ensures that n satisfies to the normalization relation 

n^nf^ = -1. 

The metric 7 induced by the spacetime metric g onto 
each hypersurface St is given by the orthogonal projector 

onto St: 

7:=g + n(8in. (2) 

Since Sj is assumed to be spacelikc, 7 is a positive defi- 
nite Riemannian metric. In the following, we call it the 
3-metric and denote by D the covariant derivative asso- 
ciated with it. The second fundamental tensor charac- 
terizing the hypersurface St is its extrinsic curvature K, 
given by the Lie derivative of 7 along the normal vector 
n: 

K := -\£r.l. (3) 

One introduces on each hypersurface St a coordi- 
nate system (x^) = {x^,x^,x^) which varies smoothly 
between neighboring hypersurfaces, so that (a;") = 



^ We use geometrized units for which G = 1 and c = 1; Greek 
^ Note that while Nakamura ct al. [30] used spherical coordinates, indices run in {0, 1, 2, 3}, whereas Latin indices run in {1, 2, 3} 

they considered Cartesian components of the tensor fields. only. 
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{t,x^,x'^,x^) constitutes a well-behaved coordinate sys- 
tem of the whole spacetime^. We denote by {d/dx°') = 
{d/dt,d/dx') = {d/dt,d/dx^,d/dx'^,d/dx^) the natu- 
ral vector basis associated with this coordinate system. 
The 3-1-1 decomposition of the basis vector d/dt defines 
the shift vector ^ of the spatial coordinates (a;*): 

d 

— =Nn + l3 with n • ^3 = 0. (4) 
at 

The metric components with respect to the coordi- 
nate system {x") are expressed in terms of the lapse func- 
tion N, the shift vector components /?* and the 3-metric 
components 'ytj according to 

g^^ dx" dx" = -N'^dt'^+jij {dx' +p'dt){dx^ +p^dt). (5) 

In the 3-1-1 formalism, the matter energy- momentum 
tensor T is decomposed as 

T = En(S'n + n(SiJ + J<S'n + S, (6) 

where the energy density E, the momentum density J 
and the strain tensor S, all of them as measured by 
the observer of 4-velocity n, are given by the follow- 
ing projections: E := T^^n^'n'^, J„ := — 7„^T^yn'', 
Sai3 '■= la^lfB^I'iJ.v- By means of the Gauss and Codazzi 
relations, the Einstein field equation is equivalent to the 
following system of equations (see e.g. Eqs. (23), (24) 
and (39) of York [3]): 

R + K'^- KijK'^ = 16wE, (7) 
DjXJ -DiK = ^TTJi, (8) 



—Kij - £0Kij = -DiDjN + N[Rij-2KikK''j 

+KKij + An {{S - E)jij - 2Sij) ] . (9) 

Equation (7) is called the Hamiltonian constraint, Eq. (8) 
the momentum constraint and Eqs. (9) the dynamical 
equations. In these equations K denotes the trace of the 
extrinsic curvature: K := K^^, S := 5'^, Rij the Ricci 
tensor associated with the 3-metric 7 and R := R^ ^ the 
corresponding scalar curvature. These equations must be 
supplemented by the kinematical relation (3) between K 
and 7: 

d 

jr,lij - £filij = -2NK,j. (10) 



B. Conformal metric 

York [38] has shown that the dynamical degrees of free- 
dom of the gravitational field are carried by the conformal 
"metric" 7 defined by 

iir-=7-'^'7ij, (11) 

where 

7:=det7y. (12) 

The quantity defined by Eq. (11) is a tensor density of 
weight —2/3, which has unit determinant and which is 
invariant in any conformal transformation of 7^ . It can 
be seen as representing the equivalence class of confor- 
mally related metrics to which the 3-metric 7 belongs. 
The conformal "metric" (11) has been used notably in 
the BSSN formulation [14, 15], along with an "associ- 
ated" covariant derivative D. However, since 7 is a ten- 
sor density and not a tensor field, there is not a unique 
covariant derivative associated with it. In particular one 
has D'y = 0, so that the covariant derivative D intro- 
duced in Sec. II A is "associated" with 7, in addition to 
D. As a consequence, some of the formulas presented in 
Refs. [14], [15] or [39] have a meaning only for Cartesian 
coordinates. 

To clarify the meaning of D and to allow for the use 
of spherical coordinates, we introduce an extra structure 
on the hypersurfaces S^, namely a metric / with the fol- 
lowing properties: (i) / has a vanishing Riemann tensor 
(flat metric), (ii) / does not vary from one hypersurface 
to the next one along the spatial coordinates lines: 

§ifij = 0, (13) 

and (iii) the asymptotic structure of the physical metric 
7 is given by /: 

7»j ~ fij spatial infinity. (14) 

This last relation expresses the asymptotic flatness of the 
hypersurfaces S^, which we assume in this article. 

The inverse metric is denoted by ^: P'^fkj = 5^ j- 
We denote by 2> the unique covariant derivative associ- 
ated with /: ffe/jj = and define 

V' ■- pVj. (15) 

Thanks to the flat metric /, we can properly define the 
conformal metric 7 as 

7y:=*"*7y or 7^^- =: ^f^^.^., (16) 
where the conformal factor is defined by 



= (7)' 



* later on we will specify the coordinates (a;') to be of spherical 
type, with = r, = 9 and = ip, but at the present stage 
we keep {x^} fully general. 



^ Note that, in general one has 7^ -y'^'^i' 
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7 and / being respectively the determinant of 7 [cf. 
Eq. (12)] and the determinant of / with respect to the 
coordinates (a;'): 

/:=det/i,-. (18) 

Being expressible as the quotient of two determinants, \E' 
is a scalar field on E^. Indeed a change of coordinates 

(.x') (x' ) induces the following changes in the deter- 
minants: 7' = (det J)^7 and /' = (det J)^/, where J 
denotes the Jacobian matrix J%, := dx^ jdx^ . It is then 
obvious that 7'//' = 7//, which shows the covariancc of 
7//. Since \E' is a scalar field, 7 defined by Eq. (16) is a 
tensor field on and not a tensor density as the quantity 
defined by Eq. (11) and considered in the BSSN formula- 
tion [1, 14, 15]. Moreover, 5* being always strictly posi- 
tive (for 7 and / are strictly positive), 7 is a Riemannian 
metric on Et . Actually it is the member of the conformal 
equivalence class of 7 which has the same determinant 
as the flat metric /: 

dot7y=/. (19) 

In this respect, our approach agrees with the point of 
view of York in Ref. [40], who prefers to introduce a spe- 
cific member of the conformal equivalence class of 7 in- 
stead of manipulating tensor densities such as (11). In 
our case, we use the extra-structure / to pick out the rep- 
resentative member of the conformal equivalence class by 
the requirement (19). 

We define the inverse conformal metric 7'^ by the re- 
quirement 

likt'=Si', (20) 
which is equivalent to 

7y=*4yj Qj. yj,^^-4~y_ ^21) 

Since 7 is a well defined metric on E(, there is a unique 
covariant derivative associated with it, which we denote 
by D: Dk')ij = 0. The covariant derivatives DT and 

X>r of any tensor field T of type on Ej are related 

by the formula 

p 

r=l 

-EA^^.T--^...,.,^, (22) 

r=l 

where A denotes the following type ( 2 ) tensor field: 

:= ^7'' + Vjju - Vn,j) . (23) 

A*^,. can also be viewed as the difference between the 



Christoffel symbols*^ of A (f'^y) and those of (r'=.^.): 

A'^.. = f^ - - f^--. (24) 

ij 13 n \ ' 

The general formula for the variation of the determi- 
nant applied to the matrix writes, once combined with 
Eq. (19), 

5\nf = 5\n^l = f^ 5%, (25) 

for any infinitesimal variation 5 which obeys Leibniz rule. 
In the special case 6 = 'Di., we deduce immediately that 

f^Vklrj = 2A' = 0. (26) 

A useful property of D is that the divergence with 
respect to it of any vector field V coincides with the 
divergence with respect to the flat covariant derivative 
■D: 

bkV^ = VkV'^ . (27) 

This follows from the standard expression of the diver- 
gence in terms of partial derivatives with respect to the 
coordinates (a;'), and from Eq. (19). 

C. Conformal decomposition 

We represent the traceless part of the extrinsic curva- 
ture by 

.= ^4 ^j^ij _ ii^7»J^ . (28) 

Again, contrary to the A^^ of the BSSN formulation [14, 
15], this quantity is a tensor field and not a tensor density. 
We introduce the following related type (2) tensor field: 

A^, := %kljiA''' = (k^j - \Kii^ , (29) 

which can be seen as A^^ with the indices lowered by 7^ j , 
instead of 7^-. Both A^^ and Aij are traceless, in the 
sense that 

jijA^^ = %A^i = and j'^A,, = 7*^i,, = 0. 

(30) 

The Ricci tensor R of the covariant derivative D (as- 
sociated with the physical 3-metric 7) is related to the 
Ricci tensor R of the covariant derivative D (associated 
with the conformal metric 7) by: 

Rij = - 2D,Dj^ + 4£),$ Dj<^ 

-2 (^D'^Dk^ + 2bk^b^^^ (31) 



Recall that, while Christoffel symbols do not constitute the com- 
ponents of any tensor field, the difference between two sets of 
them does. 
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where 



$ := In * 



(32) 



and we have introduced the notation [in the same spirit 
as in Eq. (15)] 

& := fWj. (33) 
The trace of Eq. (31) gives 

i? = *-''(i?-8i>fcI)'=$-8L>fc$ (34) 

where we have introduced the scalar curvature of the met- 
ric 



(35) 



An equivalent form of Eq. (34) is i? = '^~'^R - 

^^-^DkD^^i, which agrees with Eq. (54) of York [3]. 

Thanks to Eq. (34), the Hamiltonian constraint (7) can 
be re-written 

(36) 



This equation is equivalent to Eq. (70) of York [3]. The 
momentum constraint (8) becomes 



(37) 



which agrees with Eq. (44) of Alcubierre et al. [41] in 

the special case of Cartesian coordinates (these Authors 
are using the quantity $' = $ + 1/12 In/, with / = 1 in 
Cartesian coordinates). 

The trace of the dynamical equation (9) [combined 

with the Hamiltonian constraint (7)] gives rise to an evo- 
lution equation for the trace K of the extrinsic curvature: 



dK 
'dt 



- i3^bkK = -^-'^ {b^b^N + 2bk^b^N^ 



+N 



4n{E + S) + AkiA^^ + ^ 



whereas the traceless part of Eq. (9) becomes 



(38) 



dA^^ 



+4:{D'^D^N + D^^D'N] - i 



N 



R + sbk^b''^^ + 8bk^b''N 



+N 



r 



(39) 



where we have introduced the scalar field 
Q := ^^TV. 



(40) 



Q has the property to gather the second order derivatives 
of N and ^ in Eq. (39). Moreover, in the stationary case, 
it has no asymptotic monopolar term (decaying like 1 /r) , 
as discussed in [28] . An elliptic equation for Q is obtained 
by combining Eqs. (36) and (38): 



3 . 
4 



^^N ( 4ttS + ^AkiA''^ + ^ 
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The trace and traceless parts of the kinematical rela- 
tion (10) between K and 7 result respectively in 



— = (3^Dk^ + - {DuU" - NK) 



(42) 



and 



- £pf' - ^Dkp" = 2NA'^ . (43) 



III. EINSTEIN EQUATIONS IN TERMS OF 
THE FLAT COVARIANT DERIVATIVE 



It is worth to write the Einstein equations, not in terms 
of the conformal covariant derivative D, as done above, 

but in terms of the flat covariant derivative T>, because 
(i) numerical resolution usually proceeds through linear 
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operators expressed in terms of (and deals with non- 
linearities via iterations), and (ii) the Dirac gauge we aim 
to use is expressed in terms of D. 

A. Ricci tensor of D in terms of the flat derivatives 

of 7 

The Ricci tensor R of the covariant derivative D which 
appears in the equations of Sec. II C can be expressed in 
terms of the flat covariant derivatives of the conformal 
metric 7 as 

Rio = {Wn^j - VkVi% - VkVjjii) 

+]^Vkt' {Vaij + Vj^u - Van) 

-A^,A'.,. (44) 

This equation agrees with Eq. (2.17) of [14], provided it 
is restricted to Cartesian coordinates, for which 2?j — > di 
and JS}^^- — > ^^ij- After some manipulations, Eq. (44) 
can be written as 

Rio = —[I'^^VkVaij +^ikVjH'' +^jkVjH^ 

-A^A',., (45) 



where we have introduced the vector field 



if^:=p.f^=-7'='A'fe; (46) 



[the second equality results from Eq. (23)] . If we restrict 
ourselves to Cartesian coordinates [Vi di. A* j,; — > 
r*;.;), the quantity coincides with minus the "confor- 
mal connection functions" F* introduced by Baumgarte 
and Shapiro [15]: = —H^ . Moreover after some rear- 
rangements, the expression (45) for the Ricci tensor can 
be shown to agree with Eq. (22) of [15]. The motiva- 
tion for the writing (45) of the Ricci tensor traces back 
to Nakamura, Oohara and Kojima [30]; it consists in let- 
ting appear a Laplacian acting on [first term on the 
right-hand side of Eq. (45)] and put all the other second 
order derivatives of 7jj into derivatives of W. This is very 
similar to the decomposition of the 4-dimensional Ricci 
tensor which motivates the introduction of harmonic co- 
ordinates; note that in general the principal part of the 
Ricci tensor contains 4 terms with second-order deriva- 
tives of the metric; we have only 3 in Eq. (45) because 
det 7ij = /. 

Starting from Eq. (45), we obtain, after some compu- 
tations, an expression of the Ricci tensor in terms of the 
flat covariant derivatives of 7'-' , instead of 7jj : 



I 

+f' Im/^^fcl"" -Dnl'' + y '7fcnl?i7'"" ^mf + ^f'^ '^fe7mn ^^(7™" J • (47) 



If we restrict ourselves to Cartesian coordinates, the 
terms with second derivatives of 7*-' , i.e. the flrst three 
terms in the above equation, agree with Eq. (12) of [42]. 

The curvature scalar R defined from the Ricci tensor 
R by Eq. (35) is basically minus the flat divergence of H 
plus some quadratic terms: 

R = -VkH"^ + ^^'''Vkf^Vijij - ^^'''Vkj'^Vj^u- (48) 

B. Definition of the potentials /i'-* 

We will numerically solve not for the conformal metric 
7 but for the deviation h of the inverse conformal metric 
7'-' from the inverse flat metric, deflned by the formula 



h is a, symmetric tensor field on St of type (^) ("twice 
contravariant tensor" h^^) and we will manipulate it as 
such, without introducing any bilinear form ("twice co- 
variant tensor" hij) dual to it. 

The fiat covariant derivatives of h coincide with those 
of 7'^ : 2>fc7'-' = Vkh^^ . In particular the vector field H 
introduced in Eq. (46) is the fiat divergence of h: 

W = Vjh'^. (50) 

Thanks to the splitting (49), we can express the differ- 
ential operator j'^'-'Dk'Di which appears in the equations 
listed in Sec. Ill A as 'y'^^VkVi = A + h''''VkVi , where A 
is the Laplacian operator associated with the flat metric: 



(49) 



A := f'^'VkVi = VkV^. 



(51) 
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C. Einstein equations in terms of h and "D 



of the dynamical Einstein equations leads to 



Inserting Eq. (48) into the combination (41) of the 
Hamiltonian constraint and the trace of the spatial part 



J 



kl 



dK 
'dt 



+^2 



1 



1 



(52) 



r 



The momentum constraint (37) writes 



VjA'^+A' ^lA''^ +6A'^Vj^--f^VjK = Sw^^r, (53) 



The combination (38) of the trace of the dynamical Ein- 
stein equations with the Hamiltonian constraint equa- 
tions becomes 



with the following expression for A*^;, alternative to 
Eq. (23): 

(54) 

Taking into account property (27), the trace relation (42) 
can be expressed as 



'dt 



H^Vk^ = i {Vkl3^ - NK) . (55) 



BK 

— - p'^VkK = -*-*(A7V + h^^VkViN + H^VkN 
at ^ 



Aw{E + S) + AkiA'^' 



K 



21 



(.56) 



After some computations, the tracclcss kincmatical re- 
lation (43) and the traceless part (39) of the dynamical 
Einstein equations become respectively 



— - £ph'^ - -Vkp" h^' = 2NA'^ - {L0y\ 

fv'h^'' + V'W'^ - V^h'^ - '^H^r ) 'DkQ + s'^ 



(57) 



(58) 



where .S'-' is given by 



N 



(59) 



with 



* 2 



(60) 
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:= =-j''^Vkh" 
4 



(61) 



r 



Finally the notation {L(3y^ in Eq. (57) stands for the 
conformal Killing operator associated with the flat metric 
/ and applied to the vector field /3: 



(62) 



The writing (58) with the introduction of S"^^ by Eq. (59) 
is performed in order to single out the part which is 
linear in the first and second derivatives of h^^ (a term 
like h'^'^VkVih'^^ or h^'^V^h'-^'DiQ being considered as non- 
linear). In particular the quantities R^^ and i?* arise 
from the decomposition of the Ricci tensor (47) and Ricci 
scalar (48) in linear and quadratic parts: 

f fey'Ew = ^ {Ah'^ - V'W - V^H') + Ri^, (63) 



R 



(64) 



Conscqiicntly 5*-' contains no linear terms in the first and 
second-order spatial derivatives of h^K Regarding the 
time derivatives of h^^ (encoded in A^^), it contains only 
one linear term, in NKA^^ . Note also that the covariant 
form 7y of the conformal metric which appears in the 

expressions of lH'i and i?* is the inverse of the matrix 7*-' , 
and therefore can be expressed as a quadratic function of 
/i'-' , thanks to the fact that 7 = /. 



IV. MAXIMAL SLICING AND DIRAC GAUGE 
A. Definitions and discussion 

Let us now turn to the choice of coordinates, to fully 
specify the PDE system to be solved. First regarding the 
foliation St, we choose maximal slicing: 



(65) 



This well-known type of slicing has been introduced by 
Lichnerowicz [43] and popularized by York [3, 23]. It is 
often disregarded in 3-D numerical relativity because it 
leads to an elliptic equation for the lapse function (cf. 
discussion in Sec. I A). However it has very nice proper- 
ties: beside the well-known singularity avoidance capa- 
bility [44], it has been shown to be well adapted to the 
propagation of gravitational waves [14, 23]. 

Regarding the choice of the three coordinates (cc*) on 
each slice S^, we consider the Dirac gauge. In Dirac's 
original definition [21], it corresponds to the requirement 



d 



(V 



0. 



(66) 



This writing makes sense only with Cartesian type coor- 
dinates. In order to allow for any type of coordinates, we 
define the generalized Dirac gauge as 



1/3 



0. 



(67) 



Obviously this covariant definition is made possible 
thanks to the introduction of the flat metric / on S^. 
We recognize in Eq. (67) the flat divergence of the con- 
formal metric: 



0. 



(68) 



Since Vjf^^ = 0, this condition is equivalent to the van- 
ishing of the flat divergence of the potential /i*^ : 



Vjh'^ = 0, 



(69) 



Another equivalent definition of the Dirac gauge is re- 
quiring that the vector H vanishes [cf. Eq. (46)]: 



W = 0. 



(70) 



As discussed in Sec. I A, the Dirac gauge has been con- 
sidered as a candidate for a radiation gauge by Smarr 
and York [23] but disregarded in profit of the minimal 
distortion gauge which allows for any choice of coordi- 
nates in the initial slice. On the contrary Dirac gauge 
fully specifics (up to some boimdary conditions) the co- 
ordinates in the slices S^, including the initial one. This 
property allows the search for stationary solutions of the 
proposed system of equations. In particular this allows to 
get quasi-stationary initial conditions for the time evolu- 
tion. In this respect note that the numerous conformally 
flat initial data computed to date (see Ref. [1] for a re- 
view) automatically fulfill Dirac gauge, since the confor- 
mal flatness of the spatial metric 7 is equivalent to the 
condition h = 0. 

Another strong motivation for choosing the Dirac 
gauge is that it simplifies drastically the principal linear 
part of the Ricci tensor R. associated with the conformal 
metric: as seen on Eq. (47) or Eq. (60), this Ricci tensor, 
considered as a partial differential operator acting on h 
reduces to the elliptic term j'^^VkVih^^ in that gauge. 
Consequently, the second order part of the right hand 
side of Eq. (58) reduces to a flat Laplacian Aft*-' . This 
reduction of the Ricci tensor to a Laplacian has been 
the main motivation of the promotion of H as an inde- 
pendent variable in the BSSN formulation [14, 15]. A 
related property of the Dirac gauge is that thanks to it, 
the curvature scalar R of the conformal metric does not 
contain any second order derivative of 7'-' [set = in 
Eq. (48)]. 

Note that although Dirac gauge and minimal distor- 
tion gauge differ in the general case, both gauges result 
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asymptotically in transvcrsc-tracclcss (TT) coordinates 
(cf. Sec. IV of Ref. [23]), which are well adapted to 
the treatment of gravitational radiation. Both gauges 
are analogous to Coulomb gauge in electrodynamics. In 
1994, Nakamura [45] has used a gauge, called pseudo- 
minimal shear, which is related to the Dirac gauge, for 
it writes {d-jij / dt) = 0, while Dirac gauge implies 
Vj{d"fij/dt) =0. Note however that this pseudo-minimal 
shear does not fix the coordinates on the initial time slice, 
contrary to Dirac gauge: as the minimal distortion con- 
dition, it only rules the time evolution of the coordinate 
system. The exact Dirac gauge has been employed re- 
cently in two numerical studies, by Kawamura, Oohara 
and Nakamura [46], who call it the pseudo-minimal dis- 
tortion condition, and by Shibata, Uryu and Friedman 
[47]. 

Finally let us mention that Andcrsson and Moncrief 
[48] have shown recently that the Cauchy problem for 
3-1-1 Einstein equations is locally strongly well posed for 
a coordinate system quite similar to maximal slicing + 
Dirac gauge, namely constant mean curvature (K ~ t) 

and spatial harmonic coordinates {Vj 

0). 



B. Einstein equations within mcLximal slicing and 
Dirac gauge 

Thanks to the choices (65) and (70), the combination 
(52) of the Hamiltonian constraint and the trace of the 
spatial part of the dynamical Einstein equations simpli- 
fies somewhat 



3 , 
4 



-1-2*^ 



R 



(71) 



where we have let appear the quadratic quantity i?» de- 
fined by Eq. (61). Note that thanks to Dirac gauge, the 
right hand side of the above equation does not contain 

any second order derivative of h^^ . 

The momentum constraint (53) becomes 



Now, taking the (flat) divergence of Eq. (57) and using 
the fact that d /dt commutes with Dj, thanks to property 
(13), the Dirac gauge leads to an expression of the diver- 
gence of A^^ which does not contain any time derivative 
of the shift vector nor any second derivative of K^^ : 



A"-' 1 
^ N ' 2N 



(73) 



Inserting this relation into the momentum constraint 
equation (72) results in an elliptic equation for /3: 

A/3' + \v' {VjpJ) = IGTrTV^^J^ + 2A^^VjN 



kl 



-12NA'^Vj^ - 2NA' A 



(74) 



Thanks to maximal slicing, the kinematical trace rela- 
tion (55) reduces to 



5$ t 1 i. 

at 6 



(75) 



The combination (56) of the trace of the dynamical Ein- 
stein equations with the Hamiltonian constraint equa- 
tions becomes an elliptic equation for the lapse function: 



AN = '^^N 



4Tr{E + S) + AkiA 



kl 



- h'^^VkViN 



-2Dk^D''N. 



(76) 



In Dirac gauge + maximal slicing, the time evolution 
system (57)-(58) becomes 

— £i3h'' - -'Dkf3'' h'^ = 2NA'3 - {Lfiy^ (77) 

dt 3 2*4 ^ 

{V'h^^ + V^h'^ - V^h'^) VkQ, (78) 



VjA'^ + A' 



kl 



Aki _^ QA^JX)J<^ = 87r*V*. (72) where S'^ is slightly simphfied to 



+d,bk^b''N'j^' 



■2N 



%iA'''A'^ - An ( ^'^S'^ - -Sf^ 



(79) 
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with 



(80) 



The quadratic term i?^ in Eq. (79) is unchanged and is 
given by Eq. (61). The Lie derivatives along the shift 
vector field /3 which appear in Eqs. (77) and (78) can be 
expressed in terms of the flat covariant derivative J? by 
the standard formula: 



(81) 
(82) 



C. Wave equation for /i'-' 



As discussed in Sec. IV A, one of the main motivations 
for using Dirac gauge is that it changes the second order 

operator acting on h^^ in Eq. (78) to a mere Laplacian. It 
is therefore tempting to write the first order time evolu- 



tion system (77)- (78) as a (second order) wave equation 
for h'-^ . Note that the first order operator d/dt — £p 
which appear on the l.h.s. of the system (77)- (78) is 
nothing but the Lie derivative along the vector Nn. Its 
square is 



d 



dt "^y 5*2 

with the short-hand notation 



13' := 



a/? 

'dt 



(83) 



(84) 



Applying the operator d/dt—£j3 to Eq. (77) and inserting 
Eqs. (83) and (78) in the result leads to the wave equation 



9*2 



iV2 



1 



4 t / a 



^ - £0^ [LUr + 



r 



dt 



(85) 



Note that the left-hand side of the above equation con- 
tains all the second-order derivatives (both in time and 
space) of h'^^ , at the linear order. Actually the only 
second-order derivative of /i*^ on the right-hand side is 
the non-linear term h^^'Dj-'Dih'^^ contained in .S'-' via i?*'' 
[cf. Eqs. (79) and (80)]. 

Let us rewrite Eq. (85) as a flat-space tensorial wave 
equation: 



where □ denotes the d'Alembert operator associated with 
the flat metric / [cf. Eq. (51)]: 



5' A 



(87) 



(86) and u*^ is given by 



1 



^ - £p\ h'^ - ^Vk/3''h'^ + [L^r 



^ - £^ V„I3^ - ^(D,/3'=)2 



(88) 
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Note that we have not included into cr'-' the term^ 



d 



(89) 



which appears in the right-hand side of Eq. (85). Con- 
sequently this term appears explicitly in the right-hand 
side of Eq. (86). 

At a given time step during the evolution, a^^ is con- 
sidered as a fixed source in Eq. (86), so that the problem 
is reduced to solving a flat space wave equation. Since 
2> and □ commute (thanks to the time-independence of 
/), the source cr'-' -I- {L$y^ must be divergence-free in or- 
der for the solution /i'^ of Eq. (86) to satisfy Dirac gauge 
(69). This means that one must have 



(90) 



or, from the definition (62) of the conformal Killing op- 
erator and the vanishing of /'s Riemann tensor, 

A/3' + ip' (pj;3-') = -Vja^K (91) 

The above elliptic equation fully determines /3 (up to 
some boundary conditions), and therefore, by direct time 
integration, f3. This shows clearly that the shift vector 
propagates the Dirac spatial coordinates (x*) from one 
slice St to the next one. Hence we recover the tradi- 
tional interpretation of the shift vector. On the other 
side, (3 can be computed from the combination (74) of the 
momentum constraint and Dirac gauge condition. Both 
ways must yield the same result. However, from the nu- 
merical point of view, they may not be equivalent (due 
to numerical errors) and a strategy to compute the best 
value of (3 must be devised. 

Note that, since we reduce the time cvohition problem 
to a second-order wave equation for /i'^ , at each step, the 
extrinsic curvature term A^^ must be deduced from the 
time derivative of /i'-' and the spatial derivatives of the 
shift vector by inverting Eq. (77): 



^ ~ 2N 



(92) 



D. Transverse traceless decomposition 



The generalized Dirac gauge, expressed as Eq. (69), 
makes the potential h a transverse tensor field with re- 
spect to the metric /. However, the trace of h with 
respect to the metric /, 



docs not vanish in general, except in the linearized ap- 
proximation. Therefore h is not a transverse and trace- 
less (TT) tensor field. Since this latter property would 
considerably help the treatment of the wave equation, we 
perform a TT decomposition of h according to (see e.g. 
Sec. 7-4.2 of ADM [49]) 

h'l =: h'' + ^{h - V'V^4>) , (94) 

where is a solution of the Poisson equation 

A(j) = h (95) 

satisfying ^ = at spatial infinity. Then the trace of 

the term 1/2 [h — V^V^cf) on the right-hand side of 
Eq. (94) is equal to h. Moreover this term is divergence- 
free. We conclude that if h is transverse (Dirac gauge), 
then h defined by Eq. (94) is a TT tensor^: 



Vjh'^ = and fijh'^ = 0. 



(96) 



We then decompose Eq. (86) into a trace part, and a 
traceless one, to get 



nh 



(97) 



(98) 



where a := fijcr^^ and cr*-' is the traceless part of a^^ given 
by the decomposition analogous to (94): 

^ij =. ^ij + ^ pj _ jy^v'T) , (99) 

with AT = (J. Note that the quantity {L$y^ is trace-free 

by the very definition of operator L [Eq. (62)]. 

The search for the potentials h^^ can then proceed 
along the following steps: compute^ the; trace a of the 
effective source a^^ [Eq. (88)] and solve the Poisson equa- 
tion 



AT = cr, 



(100) 



with the boundary condition T = at spatial infinity. 
This leads to a regular solution for T because cr is a fast 
decaying source, due to the fact that Eq. (86) is the trace- 
less part, with respect to the metric -7, of the dynamical 
Einstein equations and that 7 ~ / asymptotically. The 
next step is to insert T and cr into Eq. (99) to com- 
pute (T*^ . Then one has to solve the TT wave equation 
(98) for h^^ . A resolution technique based on spheri- 
cal coordinates and spherical tensor components will be 
presented in Sec. VC. Using this technique, the resolu- 
tion of Eq. (98) is reduced to the resolution of two scalar 



h := fijh'^, 



(93) 



Eq. (89) holds thanks to the property (13). 



® If wc had removed the trace of h in the "standard" way, by 
defining h^^ := — ^hf^^ , the traceless part would not have 
been transverse. 
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d'Alembcrt equations. Then one may solve the scalar 
d'Alembert equation 

□(/) = T (101) 

for (f) and compute the trace h not by solving the 
d'Alembert equation (97) but directly as the Laplacian 
of 4> [cf. Eq. (95)]. Inserting h and (j) into Eq. (94) leads 
then to h'^^ . An alternative approach to get h will be dis- 
cussed in Sec. VD. It is algebraic [thus does not require 
to solve any d'Alembcrt equation like (97) or (101)] and 
has the advantage to enforce the condition on the deter- 
minant of 7^' [Eq. (19)]. 

V. A RESOLUTION SCHEME BASED ON 
SPHERICAL COORDINATES 



e; (g) • • • <8) e; (g) • • • <8) e-'^ (g) • • • (g) e-'" i8) are given by 

P . ... 
_|_V^pSr rpil---l.-ip 

Z-^ ik jl---jg 

r=l 

9 _ . . 

f'.. T'^-\ . (104) 

r=l 

where e-' := diag(l, 1/r, l/(rsin0)) is the change-of- 
basis matrix defined by Eq. (103), and the F*',. are the 
connection coefficients of f associated with the orthonor- 
mal frame (ej); these coefHcients all vanish, except for 



As discussed in Sec. IB, spherical coordinates have 

many advantages when treating neutron star or black 
hole spacetimes. Moreover, as we shall see below, the 
use of tensor components with respect to a spherical ba- 
sis allow to compute three of the metric components Y'' 
directly from the Dirac gauge condition (68) . In this sec- 
tion we therefore specialize the coordinates (:r' ) on each 
hypersurface to spherical ones. Moreover we expand 
all the tensor fields onto a spherical basis which is or- 
thonormal with respect to the flat metric. 

A. Spherical orthonormal btisis 

We introduce on St a coordinate system = (r, 9, ip) 
of spherical type, i.e. r spans the range [0, +oo), 9 the 
range [0,7r] (co-latitude angle), ip the range [0,27r) (az- 
imuthal angle) and the components of the flat metric / 
with respect to these coordinates are 

/i,- =diag(l, r^ r^sin^^). (102) 

The determinant / [Eq. (18)] is then / = sin^ 9. 

From the natural vector basis associated with the coor- 
dinates {r,9,(p), [d/dx^) = {d/dr,d/d9,d/dip), we con- 
struct the following vector fields: 

6^:=^, 60:=- — , e^-=^^^- (103) 
or r o9 r sin 9 Oif 

(e--) = {sr, ee,e^) forms a basis of the vector space tan- 
gent to T:f Moreover, this basis is orthonormal with 
respect to the flat metric /: f-f- = diag(l,l,l). Notice 

that we are denoting with a hat the generic indices 
associated with this basis, but we denote by r, 9, ip (with- 
out a hat) indices of specific components on this basis. 

Given a tensor field T of type the components 

of the covariant derivative in the orthonormal basis 



r^„„ — — F^ ^ — —r~^ F*" — —Vf — —r~^ 

'- 99 — '- rd — ' ) cptp — r(p ~ ' ' 

r\^^-r\^ = -{rte.ne)-\ (105) 



B. Resolution of elliptic equations 

1. Scalar Poisson equations 

We have to solve two scalar elliptic equations: the 
Hamiltonian constraint (combined with the trace of the 
dynamical Einstein equations) Eq. (71) for Q and the 
maximal slicing equation (76) for N . Both equations are 
not strictly Poisson equations since they contain Q and N 
on their right-hand side. Moreover the right-hand side of 
Eq. (71) is non-linear in Q (through <I> = (IniV — lnQ)/2). 
Therefore these equations must be solved by iterations, 
solving for a Poisson equation at each step. Since we are 
using spherical coordinates, it is natural to perform an 
expansion on spherical harmonics Yp(9, (p). The resolu- 
tion of the scalar Poisson equation is then reduced to the 
resolution of a system of second order ordinary differen- 
tial equations in r for each couple {£,m). We refer the 
reader to Ref. [19] for further details. 



2. Vector elliptic equation for the shift 

As we have seen in Sec. IV B, the Dirac gauge condi- 
tion once inserted into the momentum constraint equa- 
tion gives rise to the elliptic equation (74). Using the 
derivation formula (104) with the explicit values (105) of 
the connection coefficients, we obtain the following com- 
ponents of this equation with respect to the orthonormal 
frame (ej): 
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Qj.2 
Q^2 



2 dd'^ 
r dr 
2 9/3^ 



1 



Ae^/S-- - 2/3'- - 2 



r dr \ ^'^^ 86 

2 813^ 



r 8r 



81^ 
86 

sin^( 



2 8(3^ 



tan sin 6 8(p 
^cos6 8P^' 



159 

3 9^ 



sh^ 6 dip 
cos 6 8/3^ 



sin 6 8ip sin^ ^ 8(p sin^ ^ / 3r sin 6 8ip 



80 



(106) 
(107) 
(108) 



where Ag^ denotes the angular Laplacian: 



8^ 



+ 



8 



+ 



8^ 



86'^ tan 6 86 sm^ 6 8ip'^ 



(109) 



S{(3y stands for the right-hand side of Eq. (74) and 
:= VkP^ is the divergence of (3 with respect to the 
flat connection 27. In terms of the components with re- 
spect to the orthonormal frame (ej), it reads 



e 



813'' 2(3'' 1 /5/3' 



8r 



1 dp^ 

86 tan^ sin0 8ip 



■ (110) 



As for the scalar elliptic equations for Q and TV discussed 
above, the right-hand side S{(3Y of Eqs. (106)-(108) de- 
pend (linearly) on /3, both explicitly and via A'^ [see 
Eqs. (74) and (92)]. Thus an iterative resolution must be 
contemplated. 

Equations (106)-(108) constitute a coupled system, 
since each equation contains all the components of (3. 
However, we can decouple the system by proceedings as 
follows. First, taking the (flat) divergence of this vector 
system, and taking into account that T> and A commute 
(flat metric), we get a scalar Poisson equation for 6 only: 



AQ = -Vf^S{(3f 



(111) 



Assuming this equation is solved for 8. we use Eq. (110) 
to replace the terms containing angular components in 
Eq. (106) to get a decoupled equation for /3'': 



gj.2 



+ 



4 8P'' 2/3'- 
r 8r 



1 



A9^/3'- 



2 

3 or r 



(112) 



This equation can be solved by expanding /3'' in spherical 
harmonics. An alternative approach is to set 

X := r/3'- (113) 

which reduces Eq. (112) to an ordinary Poisson equation: 

Ax = rS{l3Y-^^+2e. (114) 

This is not surprising since x is actually a scalar field on 
^t- X = fiji'^l^^ J where r denotes the "position" vector 
field: 



(115) 



(x, y, z) and (e^, Cy, e^) being respectively the Cartesian 
coordinates and Cartesian frame canonically associated 
with the spherical coordinates (r, 9, ip). Indeed, contrary 
to Br, which is singular at the origin r = 0, r is a regular^ 
vector field [this is obvious from the second equality in 
Eq. (115)]. Being the scalar product of (3 and r (with 
respect to /), x is then a regular scalar. 

Let us now discuss the resolution of the angular part. 
We introduce a poloidal potential and a toroidal po- 
tential /i such that /3 is expanded as (see also § 13.1 of 
Ref. [55] and § A.2.a of Ref. [56]): 

(3 = 13'' Br + [rT>r] - (e^ • T>ri) r] + rx T>ijl, (116) 

where the scalar product and the vectorial product arc 
taken with respect to the flat metric /. Note that the 
terms containing r} and n are by construction tangent to 
the sphere r = const and that r x X>/i is nothing but 
the angular momentum operator of Quantum Mechanics 
applied to /i. An alternative expression is r x = 
—T) X (/ur). In term of components, Eq. (116) results in 



8rj 1 8^1 
86 sin 9 8ip 
1 8r] ^8)1 
sin 6 8(p 86 



(117) 
(118) 



It can be shown easily that the potentials rj and /U obey 
to the following relations: 



Ae^rj 
Ag^fi 



2/3"- 



9/3'- 

rB — r— — 

8r 

r ■ (I? X /3) 

(3^ 1 8(3^ 
89 tm\9 s\Ti9 8(p 



(119) 

(120) 



These formulas show that i] and fi are uniquely de- 
fined (up to the addition of some function of r). G, 
/3'' = x/r and the scalar r • (I? x /3) being expand- 
able in (scalar) spherical harmonics, Eqs. (119) and (120) 
show also that r] and /U are expandable in spherical har- 
monics Y™{6,ip). The computation of r] and ^ from 



^ As in Ref. [4], we define a regular tensor field as a tensor 
field whose components with respect to the Cartesian frame 
(e-x, ey,ez) are expandable in power series oi x, y and z. 
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the components (/?'', /?^, /J*^) can then be performed from 
Eqs. (119)-(120) by a mere division by 1) (eigen- 

value of the operator Agj^ corresponding to the eigen- 
function ¥[^{9, (p)). In the following we call this type of 
computation a quasi-algebraic one. 

By a straightforward computation, it can be shown 
that the part (107)-(108) of the original system is equiv- 
alent to the two Poisson equations 

A?? = Vs-^-i:- 121 
3 r 

An = /is, (122) 

where rjs and /is are the poloidal and toroidal potentials 
of the source S{(3) [they can thus be determined from 
S{f3) by formulas (119)-(120) with (3^ replaced by S{l3f]. 

Having reduced the complicated coupled PDE system 
(106)-(108) to Poisson type equations (111), (112), (114), 
(121) and (122), various strategies can be devised to get 
the solution. In all of them, we take advantage of the 
fact that the Poisson equation (122) for the toroidal part 
is fully decoupled from the others to solve it first and 
hence get /t. Similarly the Poisson equation (111) for the 
divergence is decoupled from the other equations. So we 



can solve it to get 9. Then we plug 8 on the right-hand 
side of Eq. (112) and solve it to get /?''. An alternative 
approach is to solve the Poisson equation (114) for x a-nd 
obtain as x/t"- Then we have the following options: (i) 
deduce rj from Eq. (119); (ii) solve the Poisson equation 
(121) to get r]. Method (ii) requires to solve an addi- 
tional Poisson equation, while method (i) requires only a 
division by —£(£+ 1) of the coefficients of spherical har- 
monics expansions, making a total of three scalar Poisson 
equations to solve the system. However method (i) in- 
volves the radial derivative of f3'^ which may result in a 
low order of differentiability of the numerical solution. 



C. Resolution of the tensor wave equation 

1. Spherical components 

By means of the derivation formula (104) with the ex- 
plicit values (105) of the connection coefficients, the ten- 
sor wave equation (98) can be written explicitly in terms 
of the components h'^^ of the TT part of h with respect 
to the orthonormal spherical basis: 





-1- 


Q^2 


2 dh"""^ 1 
r or r'^ 




9*2 
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Qj.2 


2 d\f^ 1 
r or r^ 





4/1 



re 



d0 



sin 



tan^ 



4 dh'^'f 
sin 9 dip 



86 



-2- 



89 



-2 



2h^^ + 2h'^f 
cos 9 8lf'P 



= 5"'^ (123) 
2W 



siv? 9 8ip tan 9 
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sin 9 8ip tan f 
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89 
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sin 9 8ip 
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tan2 6' 
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gj.2 
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r dr 



dip J tan 9 



AeJi'^'^ - 
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2h'^f 4 dh'^ 
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siv? 9 sin dip sin^ 9 8ip 
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(124) 



(125) 
(126) 



(127) 



(128) 
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where S^^ denotes the right-hand side of Eq. (98) : 5*-' := ct*^ + . These equations must be supplemented by 

the TT conditions [Eq. (96)], which read, in term of components with respect to (ej), 



dr 


+ 
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dJf^ 
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+ 




r 










+ 





+ 



89 



09 



1 



1 dh'''c 
sin 9 dip 

1 d¥'^ 
sin 9 dip 

1 dh.^'^ 

i 1 

sin0 dip tan^ 



hi-' 
tan^ 



+ ih^' - ft^^) 



tan^ 



= 



= 0, 



(129) 

(130) 

(131) 
(132) 



r 



As discussed in Sec. IV D, the TT conditions and the □ 
operator commute, so provided that the source S is TT, 
the solution h will also be TT. 

For the steady state case [djdt = 0) or for an implicit 
time scheme^°, we need to invert the full operator on the 
left hand side of the system (123)-(128). One immedi- 
ately notices that this system couples all the components 

A natural idea to solve the system (123)-(128) would 
be to expand h onto a a basis of tensor spherical harmon- 
ics. Notice that, contrarily to scaZar spherical harmonics, 
there are several types of tensor ones (for a review, see 
[50]). A first family has been introduced by Mathews 
[51] and Zerilli [52]; they are called pure orbital harmon- 
ics in [50] and are eigenvectors of the angular Laplacian 
(109) acting on tensors. A second family is made of pure 
spin harmonics [52, 53] which are very well suited for 
describing gravitational radiation in the radiation zone 
(where one supposes that the wave vector is parallel to 
the radial direction). However, it should be realized that 
all families of tensor spherical harmonics are based on a 
longitudinal/transverse decomposition with a notion of 
transversality different from the one used here: in our 
acceptation, transverse means divergence-free [Eqs. (69) 
and (96)], whereas in tensor spherical harmonics litera- 
ture, transverse means orthogonal with respect to the ra- 
dial vector Br- Asymptotically both notions coincide, but 
this is not the case at finite r. From the very definition 
of Dirac gauge [Eqs. (69)], it is clear that the notion of 
transversality relevant to our problem is the divergence- 
free one. As shown by Mathews [51] and explicited in the 
quadrupolar case by Teukolsky [54], it is possible a form 
linear combinations of tensor spherical harmonics which 
are divergence-free. We propose here a different route, 
which is actually simpler. We do not perform any expan- 
sion onto the tensor spherical harmonics, but use directly 
the traceless and divergence-free properties to reduce the 



With Chebyshcv spectral methods, the accumulation of col- 
location points near the boundaries implies a very severe 
Courant-Fricdrich-Levy condition and in practice prohibits ex- 
plicit schemes. 



tensor wave equation to two scalar wave equations, re- 
fiecting the two degree of freedoms of a TT symmetric 
tensor. 

Before presenting this method, let us comment upon 
another tentative of decoupling the system (123)-(128) 
that one might naively contemplate. It would consist in 
solving separately each equation (123),. ..,(128) by treat- 
ing as source the terms with h'^'- [k ^ i ov I ^ j) so that 

only an operator acting on the component Ji^^ would ap- 
pear on the left-hand side. Of course, since the other 
components of h would be present on the right-hand side, 
such a method would require some iteration. However 
this method is not applicable, due to the bad behavior 
of the truncated operator (i.e. the operator which acts 

only on hl^ in the component ij): for a regular source, it 
gives a non-regular solution. Take for instance Eq. (123) 
in the stationary case {d/dt = 0): the operator acting on 
h'''' is 



OK''' := 



2dh'' 



1 



Qj,2 

Now = x/*^^ 
scalar field on 



dr 



— {^B^W - 4/1'- 



(133) 



where x = !ik! j]h^^ r^'r^ is a regular 
see Eq. (142) below]. hJ"^ is therefore 
expandable in scalar spherical harmonics Yf"^{9,p). For 
a given {i,m), the behavior of h"' near the origin r = 
must therefore be 



(134) 



where n is some positive integer, in order for /i'"'" to be 
regular. Inserting this expression into Eq. (133) results 



Oh'' 



[n{n -l) + 2n- £{£ + 1) - 4] r"-^ yj- 



\9,ip). 
(135) 

Thus we get a regular solution of the homogeneous equa- 
tion OhJ"^ = near r = only if, for any £, there exists a 
strictly positive integer n such that -I- n — £(€-!- 1) — 4 = 
0. However in general, this last equation does not ad- 
mit any integer solution n. The generalization to the 
time-dependent case is straightforward. Moreover, even 
if r = is excluded from the computational domain (for 
example when treating black holes), a similar regularity 
problem appears in the other equations on the axis ^ = 

or TT. 
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2. Reduction to two scalar wave equations 

As mentioned above, it is possible to use the four 

TT conditions (129)-(132) to decouple the system (123)- 
(128) and to reduce it to the resolution of two scalar wave 
equations. 

A first way to proceed is to manipulate directly 
equations (123)-(132). For instance, inserting the first 
divergence-free condition (129) into (123) and using the 
traceless condition (132) results in the disappearing of 
the terms involving }f\ hT^^ , ¥^ and h'f"^: 

OT^ ar^ r (jr r'^ 

(136) 

To perform a more systematic treatment, as well as 
to gain some insight, it is worth to introduce the scalar 
product (with respect to /) of h with the position vector 
r defined by Eq. (115): 

V :=fkih"'r', (137) 

or, in term of components, 

y» = (r;^'-^r/l'■^r/^'''^). (138) 

Note that the vector field V thus defined is regular (for 
/, h and r are regular tensor fields on Ej). Prom 
the identities OV^ = fkir'ah''' + 2Vkh^^ and_ V,V^ = 
fkir^Vih''' + fijh'^ and the TT character of h, we de- 
duce immediately that the (rr, rO, rip) part of the system 
(123)-(128) with the TT conditions (129)-(132) is equiv- 
alent to the vector wave equation 

UV = fkiS'^'r^ with ViV = 0. (139) 

Let us introduce the (regular) scalar field x as the 
scalar product (with respect to /) of r and V , 

X ■■= fkir^V^ = = r^h^''. (140) 

From the identity = fkir'^nv' + 2VkV^ and the 
divergence-free character of V . we see that Eq. (139) im- 
plies the following scalar wave equation 

Hx = r'^S'"'. (141) 

Solving this equation immediately provides hJ"^ by 

if'' = ^. (142) 

Note that inserting this last relation into Eq. (136) would 
have lead directly to Eq. (141). 



wc use the same notation x as for the decomposition of the shift 
vector in Sec. V B 2, assuming that no confusion may arise. 



We then proceed as for the vector Poisson equa- 
tion treated in Sec. VB2, namely we perform the ra- 
dial/angular decomposition of V following Eq. (116)^^: 

V = Ver + [rV-q - (e^ • X>r/) r] + r x X>/x. (143) 

Combining the above equation with Eq. (138), we see 
that the potentials r} and /i are related to the components 
and /i^^ by 

" r(al"^9^) ^^^^^ 
1,^ = + (145) 

Performing the same decomposition of the source, we get: 

r \ o& smO dip ) 

= -(-^^ + ^)- (147) 
r Of of) J 

Given S^^ and S'^'^, r]§ and ii§ are computed from the 
analog of Eqs. (119)-(120) by 

= -(^ + I^-^^j-(149) 

As already discussed in Sec. VB2, the potentials rjg 
and jig are expandable in scalar spherical harmonics 
Yp{9,(p). Equations (148)-(149) are then algebraic 
{Aff^u — > —£{£ + l)u) in terms of the coefficients of the 
spherical harmonics expansion. 

The angular part of the vector wave equation (139) 
is equivalent to the following system, analogous to 
Eqs. (121)-(122) with = (since V is divergence-free) 
and V = rJi''': 

= ns-—^, (150) 

Un = ii-g. (151) 

We can see here that the equation for is fully decou- 
pled from the other equations, contrarily to that for rj 
which contains hT'^ . Actually the divergence- free con- 
dition ViV^ = relates 7] to hJ"^ by Eq. (119) (with 
y = rh''' = x/r): 



again, we use the same notation rj and /i as for the decomposition 
of /3 presented in Sec. V B 2, assuming that no confusion may 
arise. 
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This last equation can be used to compute 77, once h^^ 
has been obtained as the solution of (136) [or from the 
system (141)-(142)], instead of solving the wave equation 
(150). 

At this stage, there remains to compute the angular 
components h^^ , hP'^ and h'-^'-f . They can be deduced fully 
from the other components, by means of the TT relations 
(130)-(132). Indeed, using the traceless condition (132), 
the transverse conditions (130) and (131) can be written 
as 

-1^ + (154) 



with 



:= sin^ e ( + 3/i''' 

\ or 

:= - sin" eir^+3h''^ 
or 



Taking the angular divergence and the angular curl of 
Eqs. (153)-(154), as in Eqs. (148)-(149), we get the sys- 
tem 





rp0 


1 




(157) 


80 


-\ 

tan^ 


sm6 


' d'p 


oe 


tan^ 


1 

sin< 


? dip 


.(158) 



Again, this system is algebraic in the spherical harmonics 
representation, and therefore can be easily solved to get 
sm^eTiff and sin2 6l hP'^', after and have been 
evaluated by means of Eqs. (155)-(156). The components 
h'-^'-^ and hP'^ are then obtained by a division by sin^ 9. 
Finally is obtained by the traceless condition (132): 



scalar wave equations to be solved, Eqs. (141) and (151), 
have the form 

X^\T^{t-r,e,^) and n ^^J',.{t- r,e,ip), {WQ) 

where J^^ and are two bounded functions. Prom 
Eq. (152), we then get the following asymptotic behavior 
for the potential 77: 



-J^r,{t-r,e,ip), 



(161) 



where J-ri is a bounded function. The asymptotic behav- 
ior of the components hT"^ , Jf^ and Tf"^ follow immedi- 
ately from Eqs. (142), (144) and (145): 











(162) 


J (155) 
tan6'y* ' 




\j^i{t - r, I 




(163) 


(156) 


hTf ~ 


4j^2(i-r,i 




(164) 



where J^i and are two bounded functions. This faster 

than 0(l/r) decay shows that the (/i'''', ft''^) part 
of h does not transport any wave, as expected (cf. the 
asymptotic TT structure of Dirac gauge discussed in 
Sec. IV A). 

Thanks to the terms rdh^^/dr and rdh^'^/dr in 
Eqs. (155)-(156), it can be shown easily that the asymp- 
totic behavior of W"^ and H^f deduced from Eqs. (162)- 

(164) are 

Ji^"^ ^ --h+{t-r,9,ip) and hP'^ ^-h^{t - r,e,ip), 

(165) 

where h+ and are two bounded functions. From 
Eqs. (159), (162) and (165), one gets 



(159) 



/i^e ^ -h+(t-r,e,ip). 
r 



(166) 



In conclusion we propose to solve the tensor wave equa- 
tion (98) by solving two scalar wave equations: for x 
[Eq. (141)] and for [Eq. (151)]. h^^ is then obtained by 
dividing x by [Eq. (142)]. 77 is obtained from x by the 
quasi-algebraic equation (152). From and 77, we com- 
pute Jf^ and Iff from Eqs. (144)-(145). Then solving 
the quasi-algebraic equations (157) and (158) gives h'^'^ 
and W'f . Finally hP^ is computed by the traceless con- 
dition (159). The advantage of this procedure consists 
in solving only for two scalar wave equations which are 
linearly decoupled. This guarantees numerical stability, 
at least in the linear case. 



3. Asymptotic behavior 

Providing that the source 5^*^ is decaying sufficiently 
fast, the general asymptotic outgoing solutions of the two 



Contemplating Eqs. (165) and (166), we recover the usual 
behavior of a radiating metric in the TT gauge, h+ and 
/ix being the two gravitational wave modes. 



D. Computing the trace h by enforcing the unit 
value of the determinant of 7'^ 

Having solved the TT wave equation for h, there re- 
mains to determine the trace h = fij h^^ to reconstruct h 
by Eq. (94), and then the conformal metric 'y = f + h. 
h can be obtained by solving the scalar wave equation 
(97). However, h can also be computed in order to en- 
force a relation arising from the very definition of the 
conformal metric, namely that the determinant of the 
components 7*^ is equal to the inverse of that of the flat 
metric: detT-*-* = [cf. Eq. (19)]. It is easy to show 
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this is equivalent to the following requirement about the 
orthonormal components: 



det7*^' = 1. 

Replacing 7*^ by /'^ + h^^ , this relation writes 



1 + K"^ 



= 1. 



(167) 



(168) 



Expanding the determinant and using h = h'^'^+h^^+h'*"*' 
results in 

+ /i^^(/i'-^)2 + /i^^(/i'-^)2. (169) 

This relation shows clearly that among the six compo- 
nents h^^ only five of them are independent. The Dirac 
gauge adds three relations between the h^^ , leaving two 
independent components: the two dynamical degrees 
of freedom of the gravitational field. Equation (169) 
shows also that, at the linear order in h^^ , the condition 
detT''-' = 1 is equivalent to h = 0. 

We propose to use Eq. (169) in a numerical code to 
compute h, in order to enforce the condition (167) by 

means of the following iterative procedure: initialize h^^ 
by the TT part /i'-' obtained as a solution of the wave 
equation (98); then (i) compute h from Eq. (169); (ii) 
solve the Poisson equation (95) to get (p; (iii) insert the 
values of h and (f) into Eq. (94) to get /i'-' ; (iv) go to 
(i). In practice, this procedure converges up to machine 
accuracy (sixteen digits) within at a few iterations. 



The next equation to be solved is the TT tensor wave 
equation (98) for h, which arises from the Einstein dy- 
namical equation (78). As detailed in Sec. VC, by fully 
exploiting the TT character of h, the resolution of this 
equation is reduced to the resolution of two scalar wave 
equations for two scalar potentials x M [Eqs. (141) 
and (151)]. From x and /U one can reconstruct all the 
components of h by taking some derivatives or inverting 
some angular Laplacian (which reduces to a mere division 
by + 1) on spherical harmonics expansions). 

Then the trace h of h is determined algebraically 
through Eq. (169) which ensures that det7y = / 
[Eq. [19)]. From h and h, one reconstructs h via Eq. (94), 
at the price of solving the Poisson equation (95) for (p. 

Finally, from h, (3 and iV, one has to compute the 
conformal extrinsic curvature A^^ via Eq. (92). 

In the above scheme, the only equations which are not 
satisfied by construction are (i) Eq. (75) which relates 
the time derivative of the conformal factor ^ to the di- 
vergence of the shift vector (3 and (ii) Eq. (97) which 
is the trace part of the wave equation for h. These two 
scalar equations must however be fulfilled by the solution 
and could be used as evaluators of the numerical error. 
Alternatively, Eq. (75) could be enforced as a condition 
on VkP'' in the resolution of the elliptic equation (74) for 

In the above discussion, we have not mentioned the 
inner boundary conditions to set on some excised black 
hole. This point is discussed briefiy in Appendix A and 
will be the main subject of a future study. 



VI. FIRST RESULTS FROM A NUMERICAL 
IMPLEMENTATION 



E. A constrained scheme for Einstein equations 



A. Short description of the code 



Let us sketch the constrained scheme we propose to 
solve the full 3-D time dependent Einstein equations. 
Our aim here is not to provide; a detailed numerical al- 
gorithm, but to show how the Dirac gauge condition, in 
conjunction with the use of spherical coordinates, leads 
to a method of resolution in which the constraints are 
automatically satisfied and the time evolution equations 
are reduced to only two scalar wave equations. 

At a given time step, one has to solve the two scalar 
Poisson equations (71) and (76) to get respectively Q and 
N, and therefore the conformal factor = [Q/Ny^^. 
The Hamiltonian constraint is then automatically satis- 
fied. We have outlined the resolution technique of these 
two scalar Poisson in Sec. VBl. Let us stress here 
that a very efiicient numerical technique to solve within 
spherical coordinates scalar Poisson equations with non- 
compact support has been presented in Ref. [19]. 

Then one has to solve the vector elliptic equation (74) 
to get the shift vector /3, following the procedure pre- 
sented in Sec. VB2. The momentum constraint is then 
automatically satisfied. 



We have implemented the constrained scheme given in 
Sec. VE in a numerical code designed to evolve vacuum 
spacetimes within maximal slicing and Dirac gauge. The 
code is constructed upon the C-|— I- library Lorene [57]. 
It uses multidomain spectral methods [26, 34] to solve 
the partial-differential equations within spherical coor- 
dinates. The scalar Poisson solver is that of Ref. [19], 
whereas the vector Poisson equation for the shift is 
solved via the method (ii) presented in Sec. VB2. The 
scalar wave equations for x and /i [Eqs. (141) and (151)] 
are integrated forward in time by means of the tech- 
nique presented in Ref. [27], namely a second-order semi- 
implicit Crank-Nicholson scheme with efficient outgoing- 
wave boundary conditions. By "efficient" we mean that 
all wave modes with spherical harmonics indices £ = 0, 
1 and 2 are extracted at the outer boundary without 
any spurious reflection. This is far better than the Som- 
merfeld boundary condition commonly used in numerical 
relativity and which is valid only for the mode £ = 0. 

Various concentric shell-like domains are used, the out- 
ermost one being compactified, to bring spatial infinity 
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FIG. 1: Evolution of the /i"^"^ component of h in the plane 9 — 7r/2, between t — (upper left) and t — 8ro (lower right). The 
various snapshots are separated by a constant time interval At — ro. The size of the depicted square is 16ro, so that the wave 
extraction surface at -Rext = 8ro is given by the circle inscribed in this square. 



to the computational domain. The compactificd domain 
is employed to solve all the elliptic equations, allowing for 
the correct asymptotic flatness boundary conditions. On 
the contrary, the wave equations are solved only in the 
non-compactified domains, the outgoing-wave boundary 
condition [27] being imposed at the boundary between 
the last non-compactified shell and the compactified one. 
Further details upon the numerical code will be presented 
in a future publication. 



B. Initial data and computational setting 

We have employed the code to evolve pure 3-D grav- 
itational wave spacetimes, as in the two BSSN articles 
[14, 15]. Initial data have been obtained by means of the 
conformal thin sandwich formalism [40, 58]. The freely 
specifiable parameters of this formalism are 7, d'y/dt, K 
and dK/dt. In accordance with our choice of maximal 
slicing, we set K = Q and dK/dt = 0. Moreover, we 
use momentarily static data, d-^/dt — 0, along with a 



conformal metric 7 resulting from 

x{t^Q) = yr^exp^-^^ sin2 6isin2(p (170) 

^(^ = 0) = 0. (171) 

The constant numbers xo and ro parametrize respectively 
the amplitude and the width of the initial wave packet. 
Let us recall that, within Dirac gauge, the two scalars 
X and /i fully specify h and thus 7: (x, determine a 
unique TT tensor h according to the decomposition pre- 
sented in Sec. V C 2 and the full h is reconstructed from 
the trace h computed in order to ensure det7'^ = /~^, 
following the method given in Sec. VD. It can be shown 
that the metric defined by Eq. (170)-(171) corresponds to 
an even-parity Teukolsky wave [54] with M — 2. These 
initial data are similar to those used by Baumgarte and 
Shapiro [15] except theirs correspond to a Af = (ax- 
isymmetric) Teukolsky wave. In particular, we choose an 
amplitude xo = 10~^ similar to that in Ref. [15]. 

A total of 6 numerical domains have been used: a 
spherical nucleus of radius r = ro , surrounded by 4 spher- 
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ical shells of outer radius r = 2ro, 4ro, 6ro and 8ro, and 
an external compactified domain of inner radius r — Svq. 
The outgoing wave boundary conditions discussed above 
are set at r — 8ro, which we call the wave extraction ra- 
dius Rcxt- In particular, this means that we do not solve 
for h for r > 8ro. Consequently we set h to zero in the 
region r > SrQ. More precisely, we perform a smooth 
matching of the value of /i at r = 6ro to zero at r = 8ro. 
This means that we solving all the Einstein equations 
only for r < 6rQ. For r S [6ro,oo) we are solving the 
Einstein equations only for the lapse N, the shift vector 
/3 and the conformal factor with h set to zero in the 
r > Srg part of their source terms. We take into account 
the symmetries present in the initial data (170)-(171): 
(i) symmetry with respect to the plane 9 = Tr/2 and (ii) 
symmetry with respect to the transformation ip i—>- (p + n. 
Accordingly, the computational coordinate spans the 
interval [0, 7r/2] only and cp the interval [0,7r). In each 
domain, the following numbers of collocations points {— 
numbers of polynomials in the spectral expansions) are 
used: NrXNgxN^ = 17x9x8. The corresponding mem- 
ory requirement is 260 MB. This modest value allows the 
computation to be performed on a laptop. We have used 
two different time steps 5t = 10~^ro and 5t = 5 10~^ro, 
to investigate the effects of time discretization. 
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FIG. 3: Evolution of the ADM mass (logarithmic scale, con- 
trary to Fig. 2) for two different values of the wave extraction 
radius Tiext. 



where the integral is taken over a sphere of radius 
r = +00 and where we have adapted the original defini- 
tion [49] to general coordinates (i.e. non asymptotically 
Cartesian) by the explicit introduction of the flat metric 
/. The above integral can be re-written in terms of the 
conformal metric and conformal factor: 



3.537e-08 - 



R^^i= lOr^, dt = 0.01 r„ 
R =8r„ dt = 0.005 r„ 




FIG. 2: Evolution of the ADM mass for three different com- 
putational settings, corresponding to different values of the 
time step St and the wave extraction radius 7?oxt. 



C. Results 

The time evolution of the component h"^"^ of h is shown 
in Fig. 1. All the wave packet leaves the computational 
domain r < 8ro around t ^ 8ro and we do not notice on 
Fig. 1 any spurious reflexion. 

In order to test the code, we have monitored the ADM 
mass defined by 



A {f'lki)] dS\ (172) 



lOTT 



(173) 

Within Dirac gauge, the second term in the integrand 
vanishes identically, whereas the last one does not con- 
tribute to the integral, due to the fast decay of h (at least 
0(r~^)) implied by Eq. (169). Therefore the expression 
for the ADM mass reduces to the flux of the gradient of 
the conformal factor: 



Madm = — 



1 

2^ 



V,^dS\ 



(174) 



Hence the expression of ADM mass in Dirac gauge is 
identical to the well known expression for conformally 
flat hypersurfaces. The evolution of the ADM mass com- 
puted by means of Eq. (174) (let us recall that the sphere 
at r = oo belongs to our computational domain) is pre- 
sented in Fig. 2. For t < 3ro, one sees that the ADM 
mass is conserved, as it should be, with an accuracy of 
four digits. Moreover, Fig. 2 shows that the main source 
of error in the ADM mass is the finite value of the time 
step St. For t > 3ro, the ADM mass starts to decrease, re- 
flecting the fact that that the wave is leaving the domain 
T < ^ext- Note that by increasing the wave extraction 
radius from i?cxt = 8ro to i?cxt = lOro, we get a conser- 
vation of the ADM mass up to t ~ 5ro (dashed curved 
in Fig. 2). In Fig. 3, we present the evolution of the 
ADM mass on a longer timescale. We see clearly that, 
after remaining constant (the part shown in Fig. 2), the 
ADM mass decreases by four orders of magnitude after 
t ~ lOro (resp. t ~ 12ro) for the wave extraction radius 
Rcxt = 8ro (resp. i?cxt = lOrg). The very small value of 
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the ADM mass at late times demonstrates that all the 
wave packet has leaved the domain r < i?ext and no spu- 
rious reflection has occurred. This is due to the eflicient 
outgoing wave boundary conditions [27] set at the wave 
extraction radius. 
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FIG. 4: Relative error e [Eq. (175)] on the time derivative of 
the conformal factor in the central domain (r < ro). 
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FIG. 5: Evolution of the maximum of absolute value of the 
potential x [El- (140)] for the long term run. 

Another test is provided by Eq. (75) which relates the 
time derivative of the conformal factor ^E* to the diver- 
gence of the shift vector /3. As mentioned in Sec. VE, 
this equation must hold but is not enforced in our scheme. 
In a given numerical domain we define the relative error 
on Eq. (75) by 









max |(9$/9t| -I- max 





(175) 



where the max are taken on the considered domain. We 
represent the value of s in the domain where it is the 
largest, namely the nucleus (r < tq), in Fig. 4. We see 
that Eq. (75) is actually very well satisfied. The error is in 
fact dominated by the time discretization (second order 
scheme), and is as low as a few 10~^ for 6t = 5 10~'^ro. 



The increase of e at t ~ 4ro is spurious and is due to the 
arrival of the wave packet in the wave extraction domain 
St'o < 1" < 8ro. 

To check the long term stability of the code, we have 
let it run well after the wave packet has leaved the area 
r < 8ro, namely until t = 400ro. This very long time 
scale is similar with that used in Ref. [15] to assess the 
stability of the BSSN scheme. We found no instability to 
develop. In particular the maximum value of the poten- 
tial X remains at the round-off error value of 10~^^ that 
has been reached at t ~ 40ro (see Fig. 5). 



VII. SUMMARY AND CONCLUSIONS 

We have introduced on each hypersurface t = const of 
the 3-1-1 formalism a flat 3- metric /, in addition to the 
(physical) 3-metric 7 induced by the spacetime 4-metric 
g, in such a way that asymptotically both metrics co- 
incide. This allows us to define properly the conformal 
metric 7 and not to stick to Cartesian coordinates. A 
flat metric is introduced anyway, more or less explicitly, 
when performing numerical computations. We have writ- 
ten the 3-1-1 equations entirely in terms of the covariant 
derivative associated with the flat metric /. 

The Dirac gauge is expressed simply in terms of this 
flat metric as the vanishing of the divergence with respect 
to f of the conformal metric 7. Moreover in spherical 
components, the Dirac gauge reduces the resolution of 
the equations for 7 to two scalar wave equations. The 
remaining four components j''^ are then obtained from 
the condition det 7'^ = det /'^ and the three components 
of the Dirac condition Vj^^^ = 0. This clearly shows that 
the gravitational fleld has two degrees of freedom and this 
exhibits the TT wave behavior of the metric at infinity. 
Let us stress that the usage of spherical coordinates and 
spherical components is essential for the reduction to two 
scalar wave equations. To our knowledge, this is the first 
time that a differential gauge is used to directly com- 
pute some of the metric components, thus decreasing the 
number of PDE to be solved. Previously, this was done 
only for algebraic gauges (i.e. setting some of the metric 
components to zero). 

Contrary to e.g. the minimal distortion gauge [23] or 
the "Gamma-driver" gauge [59], the Dirac gauge com- 
pletely fix the coordinates (up to some boundary condi- 
tions) in the initial hypersurface Eq. This implies that 
initial data must be prepared within this gauge, which 
might be regarded as a drawback (for instance an an- 
alytic expression for the Kerr solution is not known in 
Dirac gauge). On the contrary, an advantage of the full 
coordinate fixing is to allow to compute stationary so- 
lutions by simply setting d/dt = in the various equa- 
tions. For instance, Shibata, Uryu and Friedman [47] 
have recently proposed to use the Dirac gauge to com- 
pute quasiequilibrium configurations of binary neutron 
stars. 

In addition to the Dirac gauge, the use of the maximal 
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slicing results in an elliptic equation for the lapse func- 
tion. Another elliptic equation for the conformal factor 
^ (or equivalently for Q := ^^N) arises from the Hamil- 
tonian constraint. The Dirac gauge itself, in conjunc- 
tion with the momentum constraint, results in an elliptic 
equation for the shift /3. The maximal slicing relates the 
divergence of /3 to the time derivative of the conformal 
factor. 

Solving the above equations implies that the four con- 
straints are fulfilled by the solution. As already men- 
tioned in the Introduction, some authors have proposed 
very recently a scheme in which the constraints, re- 
written as time evolution equations, are satisfied up to 
the time discretization errors [20]. On the contrary, in 
our scheme the constraints are fulfilled within the preci- 
sion of the space discretization errors (which can be very 
low with a modest computational cost, thanks to spectral 
methods). 

It is worth noticing that the five elliptic equations 
of the widely used Isenberg- Wilson-Mathews approxima- 
tion to General Relativity [60-62] (see also Ref. [63]) 
are naturally recovered in our scheme by simply setting 
h = 0: they arc the equations for A^, Q and /3. 

We have demonstrated the viability of the proposed 
constrained scheme by numerically computing the evolu- 
tion of a gravitational wave packet in a vacuum space- 
time. The numerical evolution has been found to be both 
very accurate and stable. We are also made confident by 
existing constrained schemes for vector equations which 
have proved to be successful: the divergence-free hydro 
scheme of Rcf. [5G] (the constraint being that the ve- 
locity field is divergence-free) and some MHD schemes 
in cylindrical coordinates [64] (the constraint being that 
the magnetic field is divergence- free) . 

In this paper we have focused on space slices with 
topology, except for Appendix A where we briefly dis- 
cuss the properties of degenerate second order operators 
and the number of boundary conditions at the surface of 
excised holes with vanishing lapse. In a future work, we 
shall focus on black hole spacetimes. 
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APPENDIX A: DEGENERATE ELLIPTIC 
OPERATORS ON A BLACK HOLE HORIZON 



coordinate vector coinciding with the Killing vector of 
stationarity) . Indeed we require arbitrary long term evo- 
lution of steady state, or quasi-steady state, black hole 
spacetimes. For classical solutions (Kerr) in usual co- 
ordinates, this requirement results in a vanishing lapse 
on the horizon (see discussion in Refs. [28, 29]). There- 
fore we excise from our computational domain a sphere 
H (or two spheres for binary systems) with A'' = as a 
boundary condition on that sphere and choose spherical 
coordinates such that r = 1 on "H 

In this case, the spatial operator acting on h in Eq. (85) 
must not be merely the Laplacian A but 

Ah'^ := NAh'^ - VkN {V'h^^ + V^h'^ - V^h'^) . (Al) 

This operator is formed by writing VkQ = "^^VkN + 
2N-^'Dk'^ in the right-hand side of Eq. (85) and gather- 
ing the VkN term with the A/i*-' one. The operator ▲ is 
degenerates because of the vanishing of N at the bound- 
ary Ti. Similarly, the operator acting on the shift vector 
/3 is degenerate on 7i (cf. Eq. (74) with A^^ given by 
Eq. (92) which contains a division by the lapse TV). Let- 
ting the unknown u be a component of /i*-' or these 
equations are of the kind 

NAu + e ViNV'u = S, (A2) 

with the associated homogeneous equation 

NAu + e V.NV'u = 0, (A3) 

where = and dN/dr > at r = 1, e = ±1 and 
S is some effective source. Since Eq. (A2) is linear, a 
solution is a linear combination of a particular solution 
and a homogeneous solution, i.e. a solution of Eq. (A3). 
In the non-degenerate case, since Eq. (A3) is of second 
order, we have two independent homogeneous solutions, 
which allow us to impose two boundary conditions. In 
the degenerate case [N = at r = 1), the number of 
regular homogeneous solutions depends upon the sign of 
e: two for e = — 1 and only one for e = To illustrate 
this, let us consider the following one-dimensional second 
order equation analogue to Eq. (A3) with x = r — 1: 

(Pu du ^ . , .„ , . 

x^ -I- e— = 0, with X G [0, 1]. (A4) 

The involved second-order operator is clearly degenerate 
at .T = 0. For e = — 1, we have two independent homoge- 
neous solutions: 

U\{x) = const and ^2(2^) — x'^, (^5) 

whereas for e = 1, the two independent homogeneous 
solutions are 

ui(a;)= const and U2{x)=\nx. (A6) 



In our view, a numerical scheme for black hole 

spacetimes should recover known stationary solutions in 
coordinate-time independent form (i.e. with the d/dt 



For a binary system, wc introduce two coordinate systems, each 
centered on one hole, cf. [25] 
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The last one is clearly not regular at x = 0, so that in 
this case, one can use only one homogeneous solution to 
satisfy a Dirichlet boundary condition. 

This behavior of the degenerate operator can also be 
understood by considering the parabolic (heat- like) equa- 
tion associated with Eq. (A3): 

— = NAu + e ViNV'u. (A7) 
ot 

The solution of the elliptic equation (A3) is the cigcnfunc- 
tion corresponding to the zero eigenvalue of the spatial 
operator acting on the right-hand side of Eq. (A7). In 
other words, the solution u wc search for is the relaxed 
solution of the heat-like equation (A7). When iV ^ 0, 
Eq. (A7) becomes an advection equation near r = 1, for 
which the number of boundary conditions at r = 1 is 



zero or one depending whether the "effective velocity" 
—eViN = —e is ingoing or outgoing at the bound- 

ary r = 1. 

For the spherical components of the shift vector, we 
have e = —1, so that a boundary condition can always 
be given at r = 1, in addition to the boundary condition 
at r = oo. Regarding the spherical components of the 
metric potential /i'-' , e = 1 for h'^'^ , which means that no 
boundary condition can be set at r = 1 in addition to 
h^^ = at r = 00. On the contrary, e = — 1 for the 
potential jj, introduced in Eqs. (144)-(145). These points 
shall be studied more in details in a future work. It is 
worth to mention that the boundary conditions for h^^ 
at r = 1 determine fully the coordinates within the Dirac 
gauge. 
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